function priceapprox=constprice5
global x r T t K % price parameters
global sigma_square_bar

x = 5; r = 0.05; T = 1; t = 0; K=4; % x has to be nonzero in order P1 != 0 
my = 1; v = 1; ay = 10;
by = v*sqrt(2*ay); rho = 0.5;


sigma_square_bar = sigmasquarebar(my,v);

P0 = conditionedprice;
[pt2P0,pt3P0] = partialP0(my,v,x,K);
[v2,v3] = v2v3(v,rho,ay,my);

P1 = (t-T)*(v2*x^2*pt2P0 + v3*x^3*pt3P0);

priceapprox = P0+P1;
end

function [P0] = conditionedprice
global x r T t K sigma_square_bar

mu = r*(T-t)-sigma_square_bar/2*(T-t);
sigma = sqrt(sigma_square_bar*(T-t));

P0 = quadgk(@(s) max(s-K,0).*lognpdf(s/x,mu,sigma)/x,-Inf, Inf,'AbsTol',1e-12,'RelTol',1e-12);%-exp(mu+sigma^2/2)^2-exp(2*mu+sigma^2)*(exp(sigma^2)-1);

P0 = exp(-r*(T-t))*P0;
end

function [ptptP0, ptptptP0]=partialP0(m,v,x,K)

f2 = @(y)  (  (m-log(y/x)).^2  +  v^2*(m-log(y/x))  - v^2) /v^4/x^2;

f3 = @(y) -(  v^2+m-log(y/x)  )  .*  (  (m-log(y/x)).^2  +  2*v^2*(m-log(y/x))  -  3*v^2  )/v^6/x^3;

ptptP0 = quadgk(@(y) max(y-K,0).*f2(y).*lognpdf(y/x,m,v)/x,-Inf, Inf,'AbsTol',1e-12,'RelTol',1e-12);
ptptptP0 = quadgk(@(y) max(y-K,0).*f3(y).*lognpdf(y/x,m,v)/x,-Inf, Inf,'AbsTol',1e-12,'RelTol',1e-12);

end



function [v2,v3]=v2v3(v,rho,ay,my)

sigma1sq = @(s) (2+sin(s)).^2;

sigma1_square_bar = sigmasquarebar(my,v);

phiprime_sigma1_Phi = @(y) 1/v^2*quadgk(@(s) (sigma1sq(s)-sigma1_square_bar).*normpdf(s,my,v),-Inf,y,...
    'AbsTol',1e-12,'RelTol',1e-12);

space = linspace(-10,10,10^4);

avarage = 0;

for ii=1:10^4
    avarage = avarage+phiprime_sigma1_Phi(space(ii));
    disp(ii);
end

avarage = avarage - phiprime_sigma1_Phi(-10)/2-phiprime_sigma1_Phi(10)/2;

avarage = avarage * 20/10^6;

v2 = v/sqrt(2*ay)*2*rho*avarage;
v3 = v2/2;
end